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LEARNING WITH CLUSTERING STRUCTURE 


VINCENT ROULET, FAJWEL FOGEL, AT,EX ANDRE D’ASPREMONT, AND FRANCIS BACH 


Abstract. We study supervised learning problems using clustering constraints to impose structure on either 
features or samples, seeking to help both prediction and interpretation. The problem of clustering features 
arises naturally in text classification for instance, to reduce dimensionality by grouping words together and 
identify synonyms. The sample clustering problem on the other hand, applies to multiclass problems where 
we are allowed to make multiple predictions and the performance of the best answer is recorded. We derive a 
unified optimization formulation highlighting the common structure of these problems and produce algorithms 
whose core iteration complexity amounts to a k-means clustering step, which can be approximated efficiently. 
We extend these results to combine sparsity and clustering constraints, and develop a new projection algorithm 
on the set of clustered sparse vectors. We prove convergence of our algorithms on random instances, based on 
a union of subspaces interpretation of the clustering structure. Finally, we test the robustness of our methods 
on artificial data sets as well as real data extracted from movie reviews. 


1. Introduction 

Adding structural information to supervised learning problems can significantly improve prediction per¬ 
formance. Sparsity for example has been proven to improve statistical and practical performance [Bach 
et al., 2012]. Here, we study clustering constraints that seek to group either features or samples, to both 
improve prediction and provide additional structural insights on the data. 

When there exists some groups of highly correlated features for instance, reducing dimensionality by 
assigning uniform weights inside each distinct group of features can be beneficial bofh in ferms of prediction 
and inferprefafion [Bondell and Reich, 2008] by significanfly reducing dimension. This often occurs in fexf 
classificafion for example, where if is nafural fo group fogefher words having fhe same meaning for a given 
fask [Dhillon el al., 2003; Jiang el al., 2011]. 

On fhe olher hand, learning a unique prediclor for all samples can be loo resfriclive. For recommendation 
systems for example, users can be parlilioned in groups, each having differenl lasles. Here, we sludy how lo 
learn a parlilion of fhe samples lhal achieves fhe besl wilhin-group prediclion [Guzman-Rivera el al., 2014; 
Zhang, 2003] 

These problems can of course be lackled by grouping synonyms or clustering samples in an unsupervised 
preconditioning step. However such partitions mighl nol be optimized or relevanl for Ihe prediction lask. 
Prior hypolheses on Ihe partition can also be added as in Lalenl Dirichlel Allocation [Blei el al., 2003] 
or Mixlure of Experls [Jordan, 1994]. We presenl here a unified framework lhal highlighls Ihe clustered 
slruclure of Ihese problems wilhoul adding prior information on Ihese clusters. While conslraining Ihe 
predictors, our framework allows Ihe use of any loss function for Ihe prediction lask. We propose several 
optimization schemes to solve Ihese problems efficienlly. 

Firsl, we formulate an explicil convex relaxation which can be solved efficienlly using Ihe conditional 
gradienl algorilhm [Frank and Wolfe, 1956; Jaggi, 2013], where Ihe core inner step amounls to solving 
a clustering problem. We Ihen sludy an approximate projected gradienl scheme similar to Ihe Iterative 
Hard Thresholding (IHT) algorilhm [Blumensath and Davies, 2009] used in compressed sensing. While 
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constraints are non-convex, projection on the feasible set reduces to a clustering subproblem akin to k-means. 
In the particular case of feature clustering for regression, the k-means steps are performed in dimension one, 
and can therefore be solved exactly by dynamic programming [Bellman, 1973; Wang and Song, 2011]. 
When a sparsity constraint is added to the feature clustering problem for regression, we develop a new 
dynamic program that gives the exact projection on the set of sparse and clustered vectors. 

We provide a theoretical convergence analysis of our projected gradient scheme generalizing the proof 
made for IHT. Although our structure is similar to sparsity, we show that imposing a clustered structure, 
while helping interpretability, does not allow us to significantly reduce the number of samples, as in the 
sparse case for example. 

Finally, we describe experiments on both synthetic and real datasets involving large corpora of text from 
movie reviews. The use of k-means steps makes our approach fast and scalable while comparing very 
favorably with standard benchmarks and providing meaningful insights on the data structure. 

2. Learning & clustering features or samples 

Given n sample points represented by the matrix X = (xi,..., XnY' G and corresponding labels 
y = (yi,..., y„), real or nominal depending on the task (classification or regression), we seek to compute 
linear predictors represented by W. Clustering features or samples is done by constraining W and our 
problems take the generic form 

minimize Loss(y, X, W) -|- R{W) 
subject to IL G W, 

in the prediction variable W, where Loss(y, X, W) is a learning loss (for simplicity, we consider only 
squared or logistic losses in what follows), R{W) is a classical regularizer and W encodes the clustering 
structure. 

The clustering constraint partitions features or samples into Q groups Qi,, Qq of size si,..., sq by 
imposing that all features or samples within a cluster Qq share a common predictor vector or coefficient 
Vq, solving the supervised learning problem. To define it algebraically we use a matrix Z that assigns the 
features or the samples to the Q groups, i.e. Ziq = 1 if feature or sample i is in group Qq and 0 otherwise. 
Denoting V = {vi,. .. ,vq), the prediction variable is decomposed as VF = ZV leading to the supervised 
learning problem with clustering constraint 

minimize Loss(y, X, IF)-|-i?(lF) 

subject to IF = ZF, Z G {0, Z1 = 1, 

in variables IF, V and Z whose dimensions depend on whether features (m = d) or samples (m = n) are 
clustered. 

Although this formulation is non-convex, we observe that the core non-convexity emerges from a clus¬ 
tering problem on the predictors W, which we can deal with using k-means approximations, as detailed in 
Section 3. We now present in more details two key applications of our formulation: dimensionality reduction 
by clustering features and learning experts by grouping samples. We only detail regression formulations, 
extensions for classification are given in the Appendix 7.1 . Our framework also applies to clustered mul¬ 
titask as a regularization hypothesis, and we refer the reader to the Appendix 7.2 for more details on this 
formulation. 

2.1. Dimensionality reduction: clustering features. Given a prediction task, we want to reduce dimen¬ 
sionality by grouping together features which have a similar influence on the output [Bondell and Reich, 
2008], e.g. synonyms in a text classification problem. The predictor variable W is here reduced to a single 
vector, whose coefficients take only a limited number of values. In practice, this amounts to a quantization 
of the classifier vector, supervised by a learning loss. 
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Our objective is to form Q groups of features Qi,..., Qq, assigning a unique weight Vq to all features in 
group Qq. In other words, we search a predictor m G such that Wj = Vq for all 3 e Gq- This problem can 
be written 

minimize ^ YIi=i loss (jji.uFxi) + f ||m||| 
subject to w = Zv, Z G {0, ^ Z1 = 1, 

in the variables w £ W^, V £ and Z. In what follows, loss(yj, nFxi) will be a squared or logistic loss 
that measures the quality of prediction for each sample. Regularization can either be seen as a standard I 2 
regularization on w with R{w) = | ||u)|||, or a weighted regularization on v, R{v) = | Sgll^glli- 

Note that fused lasso [Tibshirani et ah, 2005] in dimension one solves a similar problem that also quan¬ 
tizes the regression vector using an li penalty on coefficient differences. The crucial difference with our 
setting is that fused lasso assumes that the variables are ordered and minimizes the total variation of the 
coefficient vector. Here we do not make any ordering assumption on the regression vector. 

2.2. Learning experts: clustering samples. Mixture of experts [Jordan, 1994] is a standard model for 
prediction that seeks to learn Q predictors called “experts”, each predicting labels for a different group of 
samples. For a new sample x the prediction is then given by a weighted sum of the predictions of all experts 
y — Pq'^q The weights pq are given by a prior probability depending on x. Here, we study a slightly 
different setting where we also learn Q experts, but assignments to groups are only extracted from the labels 
y and not based on the feature variables x as illustrated by the graphical model in Figure 1 . 



Figure 1 . Learning multiple diverse experts {left), mixture of experts model (right). The 
assignment matrix Z gives the assignment to groups, grey variables are observed, arrows 
represent dependance of variables. 


This means that while we learn several experts (classifier vectors), the information contained in the fea¬ 
tures X is not sufficient to select the best experts. Given a new point x we can only give Q diverse answers 
or an approximate weighted prediction y = Yl^=i ^'^'q^- algorithm will thus return several answers 
and minimizes the loss of the best of these answers. This setting was already studied by Zhang [2003] for 
general predictors, it is also related to subspace clustering [Elhamifar and Vidal, 2009], however here we 
aheady know in which dimension the data points lie. 

Given a prediction task, our objective is to find Q groups Qi,...,Qq of sample points to maximize 
within-group prediction performance. Within each group Qq, samples are predicted using a common linear 
predictor Vq. Our problem can be written 


1 


minimize — 

n 


^^loss {yi,v'^Xi) +^'^Sq\\Vq\\l 
q=li&Gq q=t 


( 3 ) 


in the variables V = (vi,. .., vq) £ R'^^^ and Q = (Qi, . .., Qq) such that ^ is a partition of the n 
samples. As in the problem of clustering features above, loss(yj, measures the quality of prediction 
for each sample and R{V) = | 'Sglkglli ^ weighted regularization. Using an assignment matrix 

Z £ {0,and an auxiliary variable W = (wi, ..., Wn) £ such that W = VZ"^, which means 
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Wi = Vq\l i ^ Qq, problem (3) can be rewritten 

minimize ^ Ya=i loss (vi, wf xA + f Ya=i 11^*Hi 
subject to W'^ = ZV'^, Z G {0, Ij^xQ, = 1, 

in the variables W G y g j^cixQ and Z. Once again, our problem fits in the general formulation 

given in (1) and in the sections that follows, we describe several algorithms to solve this problem efficiently. 

3. Approximation algorithms 

We now present optimization strategies to solve learning problems with clustering constraints. We begin 
by simple greedy procedures and a more refined convex relaxafion solved using approximafe condifional 
gradienf. We will show fhaf fhis laffer relaxafion is exacf in fhe case of feafure clusfering because fhe inner 
one dimensional clusfering problem can be solved exacfly by dynamic programming. 

3.1. Greedy algorithms. For bofh clusfering problems discussed above, greedy algorifhms can be derived 
fo handle fhe clustering objective. A slraighlforward sfrafegy fo group fealures is fo firsl frain predicfors as 
in a classical supervised learning problem, and fhen cluster weighfs fogefher using k-means. In fhe same 
spirif, when clusfering sample poinfs, one can alfernafe minimizafion on fhe predicfors of each group and 
assignmenf of each poinf fo fhe group where ifs loss is smallesf. These mefhods are fasl buf unsfable and 
highly dependenf on inifializafion. However, alternating minimization can be used fo refine fhe solufion of 
fhe more robusf algorifhms proposed below. 

3.2. Convex relaxation using conditional gradient algorithm. Another approach is to relax the problem 
by considering the convex hull of the feasible set and use the conditional gradient method (a.k.a. Frank- 
Wolfe, [Frank and Wolfe, 1956; Jaggi, 2013]) on the relaxed convex problem. Provided that an affine 
minimization oracle can be computed efficiently, the key benefit of using this method when minimizing a 
convex objective over a non-convex set is that it automatically solves a convex relaxation, i.e. minimizes 
the convex objective over the convex hull of the feasible set, without ever requiring this convex hull to be 
formed explicitly. 

In our case, the convex hull of the set {W : W = ZC, Z G {0, l}"^xQ^ = 1} is the entire space 

so the relaxed problem loses the initial clustering structure. However in the special case of a squared loss, 
i.e. loss{y,y) = ^{y — y)^, minimization in V can be performed analytically and our problem reduces 
to a clustering problem for which this strategy is relevant. We illustrate this simplification in the case of 
clustering features for a regression task, detailed computations and explicit procedures for other settings are 
given in Appendix 7.3. 

Replacing w = Zv in (2), the objective function in problem (2) becomes 

1 A 

(l){v,z) = —J2{yi-iZvfxiy + -\\Zv\\l 

i=l 

= —v^Z^X^XZv + -v'^Z^Zv - -y'^XZv + —y^y. 

2n 2 n 2n ^ 

Minimizing in v and using the Sherman-Woodbury-Morrison formula we then get 

min^{v,Z) = —y'^(l-XZ{Z^X^XZ + XnZ^Z)-^Z^X^)y 
v2n 

and the resulting clustering problem is then formulated in terms of the normalized equivalence matrix 

M = Z{Z'^Z)-^Z'^ 

such that Mij = 1/sq if item i and j are in the same group Qq and 0 otherwise. 
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Writing A4 = {M : M = Z{Z'^ Z) ^Z"^, Z G {0, Z1 = 1} the set of equivalence matrices for 

partitions into at most Q groups, our partitioning problem can be written 

minimize (I + MX'^'j ^ y 

subject to M £ M. 

in the matrix variable M £ Sn- We now relax this last problem by solving it (implicitly) over the convex 
hull of the set of equivalence matrices using the conditional gradient method. Its generic form is described 
in Algorithm (1), where the scalar product is the canonical one on matrices, i.e. {A, B) = Tr(A^i3). At 
each iteration, the algorithm requires solving an linear minimization oracle over the feasible set. This gives 
the direction for the next step and an estimated gap to the optimum which is used as stopping criterion. 


Algorithm 1 Conditional gradient algorithm 

Initialize Mq £ M. 
for f = 0,..., T do 

Solve linear minimization oracle 

At = argmin {N,Xil}{Mt)) (4) 

Afehull(A4) 


if gap(Mt,M*) < e then 
return Mt 
else 

Set Mt+i = Mt + at (At — Mt) 

end if 
end for 


The estimated gap is given by the linear oracle as 

gap(Mt,M,) A -(At - Mt, VV^(Mt)). 

By definition of the oracle and convexity of the objective function, we have 

-(At - Mt,X^P{Mt)) > -{M, - MuXi^{Mt)) > i^{Mt) - ij{M, 


Crucially here, the linear minimization oracle in (4) is equivalent to a projection step. This projection step is 
itself equivalent to a k-means clustering problem which can be solved exactly in the feature clustering case 
and well approximated in the other scenarios detailed in the appendix. For a fixed matrix M £ hull(Af), 
we have that 

p A -vV’(M) = ^x'^(j + Vy^(i + \xMx'^)-^x 

Zn^A nX nX 

1 

is positive semidefinite (this is the case for all the settings considered in this paper). Writing Pa its matrix 
square root we get 


argmin {N,X'ip{M)) 
Arehull(At) 


argmin Tr {N'^Vijj{M)) 

N&M 

, 1 

argmin — Tr(AP 2 Pa ) 

ngM 

argminTr((I — AjPaPa )) 

NeM 

II i 1 Il2 

argmin yPa —NP'^Wp 
N&M 

argmin min ||P 2 — ZV\\‘^p, 
z ^ 
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because N is an orthonormal projection (iV^ = N, = N) and so is (/ — N). Given a matrix W, we 
also have 

Q 

argmin ||VK — ZV\\j^ = argmin Uruj — Vq\\ 2 , (5) 

where the minimum is taken over centroids Vq and partition {Qi,... ,Qq). This means that computing 
the linear minimization oracle on VV'(M) is equivalent to solving a k-means clustering problem on 
This k-means problem can itself be solved approximately using the k-means-i-i- algorithm which performs 
alternate minimization on the assignments and the centroids after an appropriate random initialization. Al¬ 
though this is a non-convex subproblem, k-means-i-i- guarantees a constant approximation ratio on its solu¬ 
tion [Arthur and Vassilvitskii, 2007]. We write k-means(17, Q) the approximate solution of the projection. 
Overall, this means that the linear minimization oracle (4) can therefore be computed approximately. More¬ 
over, in the particular case of grouping features for regression, the k-means subproblem is one-dimensional 
and can be solved exactly using dynamic programming [Bellman, 1973; Wang and Song, 2011] so that 
convergence of the algorithm is ensured. 

The complete method is described as Algorithm 2 where we use the classical stepsize for conditional 
gradient at = A feasible solution for the original non-convex problem is computed from the solution 
of the relaxed problem using Frank-Wolfe rounding, i.e. output the last linear oracle. 


Algorithm 2 Conditional gradient on the equivalence matrix 

Input: X,y,Q,e 
Initialize Mq G M 
for f = 0,..., T do 

Compute the matrix square root P 2 of — VV’(Mo) 

Get oracle A* = k-means(P 2 ^ Q) 
if - Tr{At - MtfVil}{Mt) < e then 
return Mt 
else 

Set Mt+i = Mt + at{At — Mt) 

end if 
end for 

Z* is given by the last k-means 

V* is given by the analytic solution of the minimization for Z* fixed 

Output: V *, Z* 


3.3. Complexity. The core complexity of Algorithm 2 is concentrated in the inner k-means subproblem, 
which standard alternating minimization approximates at cost 0{tKQp), where tK is the number of alter¬ 
nating steps, Q is the number of clusters, and p is the product of the dimensions of V. However, computation 
of the gradient requires to invert matrices and to compute a matrix square root of the gradient at each iter¬ 
ation, which can slow down computations for large datasets. The choice of the number of clusters can be 
done given an a priori on the problem (e.g. knowing the number of hidden groups in the sample points), or 
cross-validation, idem for the other regularization parameters. 

4. Projected Gradient algorithm 

In practice, convergence of the conditional gradient method detailed above can be quite slow and we also 
study a projected gradient algorithm to tackle the generic problem in (1). Although simple and non-convex 
in general, this strategy used in the context of sparsity can produce scalable and convergent algorithms in 
certain scenarios, as we will see below. 
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4.1. Projected gradient. We can exploit the fact that projecting a matrix W on the feasible set 

{W:W = ZV,Ze {0, l}™xQ, Z1 = 1} 


is equivalent to a clustering problem, with 


Q 

argmin || VF — ZV\\\ = argmin Uru, — Vq\\\, 

where the minimum is taken over centroids Vq and partition {Gi, ■ ■., Gq)- The k-means problem can be 
solved approximately with the k-means++ algorithm as mentioned in Section 3.2. We will analyze this 
algorithm for clustering features for regression in which the projection can be found exactly. Writing k- 
means(l/, Q) the approximate solution of the projection, cj) the objective function and at the stepsize, the 
full method is summarized as Algorithm 3 and its implementation is detailed in Section 4.3. 


Algorithm 3 Proj. Gradient Descent 

Input: X,y,Q,e 
Initialize Wq = 0 
while \<PiWt) - G{Wt-i)\ > e do 

1 =Wt- at(V Loss(y, X, Wt) + VR{Wt)) 
[Zt+i, Vt+i] = k-means(iyj_^ 1, Q) 

Wt+i = Zt+iVt+i 

end while 

Z* and V* are given through k-means 
Output: W*,Z*,V* 


4.2. Convergence. We now analyze the convergence of the projected gradient algorithm with a constant 
stepsize at = 1, applied to the feature clustering problem for regression. We focus on a problem with 
squared loss without regularization term, which reads 

minimize ^\\Xw — y \\2 

subject to w = Zv, Z G {0,Z1 = 1 

in the variables w G v G and Z. We assume that the regression values y are generated by a linear 
model whose coefficients w* satisfy the constraints above, up to additive noise, with 

y = Xw* + 7] 

where y ~ A/’(0, u^). Hence we study convergence of our algorithm to w*, i.e. to the partition G* of its 
coefficients and its Q values. 

We will exploit the fact that each partition G defines a subspace of vecfors w, so the feasible set can be 
written as a union of subspaces. Let ^ be a partition and define 

Ug = {w w = Zv, Z G Z(G)}, 

where Z(G) is fhe sef of assignment matrices corresponding to G- Since permuting the columns of Z 
together with the coefficients of v has no impact on w, the matrices in Z(G) are identical up to a permutation 
of their columns. So, for Z G Z(G), Z(G) = {ZII, H permutation matrix}, therefore L/g is a subspace and 
the corresponding assignment matrices are its different basis. 

To a feasible vector w, we associate the partition G of its values that has the least number of groups. 
This partition and its corresponding subspace are uniquely defined and, denoting V the set of partitions in 
at most Q clusters, our problem (6) can thus be written 

minimize ^\\Xw — y\\l 
subject to we [Jg^'pUg . 
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where the variable w belongs to a union of subspaces Ug. 

We will write the projected gradient algorithm for (6) as a fixed point algorithm whose contraction factor 
depends on the singular values of the design matrix X on collections of subspaces generated by the parti¬ 
tions Q. We only need to consider largest subspaces in terms of inclusion order, which are the ones generated 
by the partitions into exactly Q groups. Denoting Vq this set of partitions, the collections of subspaces are 
defined as 

£l = {Kg, Q G Vq}, 

£2 = +^^2) (^1) ^ 2 ) G 

£3 = {UQ^VUg^VUg^, {Q\,Q‘ 2 ,,Q 3 ) ^Vq}. 

Our main convergence resulf follows. Provided fhaf fhe confracfion factor is sufficienl, if sfafes fhe con¬ 
vergence of fhe projecfed gradienf scheme fo fhe original vecfor up fo a consfanf error of fhe order of fhe 
noise. 


Proposition 4.1. Given that projection on [_}g^-phlg is well defined, the projected gradient algorithm ap¬ 
plied to (3) converges to the original w* as 


where 


I * II ^ til 4^ 

IW — Wt\\2 < P \\W 


, 1 - P* , 

2 + r,- 

1 - p 


p = 2max||/- 
lA^Sz 

V = - m.Siy.\\XIlu \\2 

71 hi^S 2 


and Hu is any orthonormal basis of the subspace lA. 


Proof. To describe fhe algorifhm we define Qt and as fhe partitions associafed respectively wifh wt 
and w* confaining fhe leasf number of groups and 

^ Wt+i /2 = Wt-X'Loss{X,y,wt) = Wt - £x'^X{wt - w*)£x'^rj 
wt+i = argmin^gy^^^^g ||u> - 

£lt = Ug, 

£it,* = lApt 

lAt,t+i,* = • 

Orfhonormal projections on Ut, Ut,* and Ut,t+i,* are given respectively by Pt, Pt^^, Pt^t+i,*- Therefore by 
definition wt G Ut, {wt, w*) G Ut,* and {wt, wt+i-,w*) G Ut,t+i,*- 
We can now confrol convergence, wifh 

||tu* - t(;t+i||2 = - u;t+i)||2 

< \\Pt+l,*{w* - Wt+i/ 2)\\2 + \\Pt+l,*{Wt+l /2 - Wt+l)\\ 2 - 

In fhe second term, as w* G [Jg^pUg and wt+i = argmin^gy^^^^^ Hm — Wtpi/ 2 \\ 2 , we have 

\\wt+l - Wt+i/ 2 \\l < Ih* - Wt+i/ 2 \\l 

which is equivalenf fo 

\\Pt-Vl,*iwt-Vl — ^^ 4 + 1 / 2)112 + IK-^ “ Pt-Vl,*)'tft+l/2\\2 ^ WPt-Vlpiu!* — Wt+l/2)\\^ + ll(-^ “ )'Zi't+l/21|2 

and fhis lasf sfafemenf implies 

\\Pt+l,*{wt+l - Wt+i/2)\\2 < \\Pt+l,*{w* - Wt+i/2)\\2- 


( 6 ) 




This means that we get from ( 6 ) 

llm* - mt+i ||2 < 2 

= 2 

< 2 

= 2 


< 2 


Pt+l,^{w* - Wt+i/2)\\2 


Now, assuming 


Pt+i^^{w* -Wt- -X'^X{w* - Wt) - -X'^r ])\\2 
n n 

Pt+l4l - -X^X){w* - Wt )\\2 + -||Pm,*(^^^)l |2 
n n 

Pt+i4l - -X^X)Pt4w* - wt )\\2 + -||Pt+i,*(X'^r ?)||2 

n n 

Pt+i,,{I - -X^X)PMw* - Wth + -\\Pt+i,.X^h\\rih. 

n n 


2\\Pt+i,,{I - -X^X)Pt,42 < p 

n 




n 


\2 ^ 


< 12 


(V) 

( 8 ) 


and summing the latter inequality over t, using that mo = 0 , we get 

1 — n* 

Ik* - Wth < P*lk*ll2 + q- —hlvh- 

1 - p 

We bound p and v using the information of X on all possible subspaces of <^2 or £^ 3 . For a subspace U ^ 62 
or £^ 3 , we define Pu the orthonormal projection on it and Hu any orthonormal basis of it. For ( 8 ) we get 

||Pt+i,*X ^||2 = ||XPt+i ,*||2 < max ||XPiy ||2 = max ||Xnw|| 2 , 

which is independent of the choice of Hu- 

For (7), using thsitUt,* C Ut,t+i,* and£ft+i,* C Upt+i,*, we have 

||Pi+l,*(/-X^X)Pt,*||2 < ||Pt,t+l,*(/- ^X^X)Pt,i+i,,||2 

n 

< m?t^\\Pu{I-^X^X)Pu \\2 
u^E-i n 

= max ||nw(/ - ^Uj;X^XUu)Uj ;\\2 
= max ||/-117^X^X11^/112, 

which is independent of the choice of 11^ and yields the desired result. ■ We now show that p and z/ derive 


from bounds on the singular values of X on the collections S 2 and E^. Denoting Smm(^) and Smax(^) 
respectively the smallest and largest singular values of a matrix A, we have 

max ||Xn //||2 = maxsmax(-^n//), 
hi^S 2 hi^S 2 


and assuming G 3 and that 

1 ^ E: Smin 


(XUu 

V 


< 5 


_£ ^max 


(XHu\ 

V xk ) 


< 1 + (5, 


for some <5 > 0, then [Vershynin, 2010, Lemma 5.38] shows 


1 . 


|/- -n^X^Xnw ||2 < 3max{<5,k}- 
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We now show that for isotropic independent sub-Gaussian data Xi, these singular values depend on the 
number of subspaces of <5i, N, their dimension D and the number of samples n. This proposition reformu¬ 
lates results of Vershynin [2010] to exploit the union of subspace structure. 

Proposition 4.2. Let £ 1 , 82 ,£z be the finite collections of subspaces defined above, let D = dim(C/) 

and N = Card(iSi). Assuming that the rows Xi of the design matrix are n isotropic independent sub- 
gaussian, we have 

^ max 11 X 11^^112 < 1 + (52 + e and max ||/--nJx'^Xniy ||2 < SmaxjJs + e, (^3 + e)^}, 
yjn U&S 2 UeSs n 

with probability larger than 1 — exp(—ce^n), where dp = any orthonormal 

basis oflA and Cq, c depend only on the sub-gaussian norm of the Xi. 

Proof. Let us fix ^ G £p, with p = 2 or 3 and one of its orthonormal basis. By definition of £p, 
AmiifA) < pD. The rows of XHu are orthogonal projections of the rows of X onto U, so they are still 
independent sub-gaussian isotropic random vectors. We can therefore apply [Vershynin, 2010, Theorem 
5.39] on Xliu G s > 0, with probability at least 1 — 2 exp(—cs^), the smallest 

and largest singular values of the rescaled matrix written respectively and Smax () are 

bounded by 


1 -Co 



n 


fXUu\ ^ 
^ ^min 1 !— I ^ 





pD s 

H-^5 

n \/n 


( 9 ) 


I /— 1-2: ^max 

\ yjTi J 

where c and Cq depend only on the sub-gaussian norm of the Xj. Now taking the union bound on all subsets 
of £p, (9) holds for any U ^ £p with probability 


1 - 2 


exp(—cs^) > 1 


_,(e_NV 
P J 


exp(—cs^) 


> 1 — 2 exp(l -|- p log(X) — cs^) 

Taking s = _|_ eyfi, we get for all U G £p, 

(XJ\u\ 


1 — 5p — e < s 


^min 


< S 


—2 ^max 


\ s/n ) 

with probability at least 1 — exp(—ce^n), where 5p = Cq- '^ 




+ 


l-l-p logfA'') 


. Therefore 


max llXfli^lb < 1 + (52 + e. 


Then [Vershynin, 2010, Theorem 5.39] yields 


1 . 


max ||/- XYiu \\2 < 3max{(53 -f e, ((53 -f e) }, 

Tl 


hence the desired result. 


Overall here. Proposition 4.1 shows that the projected gradient method converges when the contraction 
factor p is strictly less than one. When observations x* are isotropic independent sub-gaussian, this means 


Co\ — < - and 
V n 3 


l-|-31og(X) 1 

cn 3 


which is also 


n = £l{D) and n = 0(log(X)) 
10 


( 10 ) 

















The first condition in (10) means that subspaces must be low-dimensional, in our case D = 2,Q and we 
naturally want the number of groups Q to be small. The second condition in (10) means that the structure 
(clustering here) is restrictive enough, i.e. that the number of possible configurations, N, is small enough. 

As we show below, in the simple clustering case however, this number of subspaces is quite large, growing 
essentially as 

Proposition 4.3. The number of subspaces N in £i is lower bounded by 

N > 

Proof. El is indexed by the number of partitions in exactly Q clusters, /.e.the Stirling number of second 
kind Iq}. Standard bounds on the Stirling number of the second kind give 

\{Q^ + Q + - 1 < 1^1 < \{ed/Q)^Q^-^. (11) 

hence N > ■ 

This last proposition means that although the intrinsic dimension of our variables is of order D = 2>Q, the 
number of subspaces N is such that we need roughly n > 3dlog{Q), i.e. approximately as many samples 
as features, so the clustering structure is not specific enough to reduce the number of samples required by 
our algorithm to converge. On the other hand, given this many samples, the algorithm provably converges 
to a clustered output, which helps interpretation. 

As a comparison, classical sparse recovery problems have the same structure [Rao et ah, 2012], as k- 
sparse vectors for instance can be described as {ru : tu = Zv, Z'^1 = 1} and so are part of a “union 
of subspaces”. However in the case of sparse vectors the number of subspaces grows as d^ which means 
recovery requires much less samples than features. 

4.3. Implementation and complexity. In our implementation we use a backtracking line search on the 
stepsize at that guarantees decreasing of the objective. At each iteration if 

Wt+i = k-means {Wt - at(V Loss(y, X, Wt) + VR{Wt)), Q) 

decreases the objective value we take Wt+i = Wt+i and we increase the stepsize by a constant factor 
at+i = aat with a > 1. If Wt+i increases the objective value we decrease the stepsize by a constant factor 
at ■<— bat, with h < 1, output a new Wt+i and iterate this scheme until Wt+i decreases the objective value 
or the stepsize reaches a stopping value e. We observed better results with this line search than with constant 
stepsize, in particular when the number of samples for clustering features is small. 

Complexity of Algorithm 3 is measured by the cost of the projection and the number of iterations until 
convergence. If approximated by k-means-i-i- the projection step costs 0{tKQp), where tx is the number 
of alternating steps, Q is the number of clusters, and p is the product of the dimensions of V. When 
clustering features for regression, the dynamic program of Zhang [2003] solving exactly the projection step 
is in 0{d‘^Q) and ours for A:-sparse vectors, detailed in Section 5.1, is in 0{k^Q). We observed convergence 
of the projected gradient algorithm in less than 100 iterations which makes it highly scalable. As for the 
convex relaxation the choice of the number of clusters is done given an a priori on the problem. 

5. Sparse and clustered linear models 

Algorithm 3 can also be applied when a sparsity constraint is added to the linear model, provided that the 
projection is still defined. This scenario arises for instance in text prediction when we want both to select a 
few relevant words and to group them to reduce dimensionality. Formally the problem of clustering features 
(2) becomes then 

minimize ^ Ya=i loss [yi^vFxf) + ^||tu||i 

subject to u; = SZv, S G {0, Z G {0, Ij^xQ, S'^1 = 1,Z1 = 1, 
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in the variables w G v G S and Z, where S' is a matrix of k canonical vectors which assigns 
nonzero coefficients and Z an assignment matrix of k variables in Q clusters. 

We develop a new dynamic program to get the projection on /c-sparse vectors whose non-zero coefficients 
are clustered in Q groups and apply our previous theoretical analysis to prove convergence of the projected 
gradient scheme on random instances. 


5.1. Projection on fe-sparse Q-clustered vectors. Let W be the set of A:-sparse vectors whose non-zero 
values can be partitioned in at most Q groups. Given x G we are interested in its projection on W, 
which we formulate as a partitioning problem. For a feasible w G W, with Go = {i : Wi = 0} and 
Gi, ■ ■ ■, Gq', with Q' < Q, the partition of its non-zero values such that Wi = Vq if and only if i G Gq, the 
distance between x and w is given by 


X-w\\l = 

i&Go <?=! i&Gq 


The projection is solution of 

minimize ZieGo + S?=i .^2) 

subject to Card ^U^=i ^ k, 0 < Q' < Q, 

in the number of groups Q', the partition G = (Go, ■ ■ ■, Gq') of { 1 ,..., d} and v G . For a fixed number 
of non-zero values k', fhe objective is clearly decreasing in fhe number of groups Q', which measures fhe 
degrees of freedom of fhe projection, however if cannof exceed k'. We will use Ibis argumenf below fo gef 
fhe besf parameter Q'. For a fixed parfifion G, minimization in v gives fhe barycenfers of fhe Q' groups, 
k-q — 't' '^i&Gq Irissrting fhem in ( 12 ), fhe objective can be developed as 


i^Go <?=! i&Gq *=1 <?=1 


Splitting Ibis objecfive befween posifive and negafive barycenters, we gel lhal fhe minimizer of (12) solves 

maximize ^ ^ Sqfj^l 

subjecf fo Card < k, 0 < Q' < Q, ^ 

in fhe number of groups Q' and fhe partition G = (Go, ■ ■ ■, Gq') of { 1 ,..., d}, where p-q = Xi- 

We fackle Ibis problem by finding fhe besf balance befween fhe Iwo lerms of fhe objecfive. We define 
/_ (j, q) fhe oplimal value of ^ Spfj.p when picking j poinfs clusfered in q groups forming only negafive 
barycenfers, i.e. fhe solution of fhe problem 

maximize ^pkl 

subjecf fo /Xp = ^ EisPp 3 :* < 0 ( 2 - g)) 

Card ^Up=i Gp) = j, 

in fhe partition G = {Go, ■ ■ ■, Gq{ of d}. We define /+(y, q) similarly excepf lhal if conslrainls 

barycenfers fo be posifive. Using remark above on fhe paramefer Q’ , problem (13) is Ihen equivalenl fo 


maximize /_ (j, q) + f+ (k' - j, Q' - q) 

subjecf fo 0<k'<k, 0<j< k', Q' = min(A:', Q), 0 < g < Q\ 


in variables j, k' and q. 

Now we show lhal /_ and /+ can be computed by dynamic programming, we begin wilh /_. Remark 
lhal (P-{j, q)) is a partitioning problem on fhe j smallesl values of x. To see Ihis, lei S- C {1,..., d} be 
fhe oplimal subsel of indexes faken for iP-(j, q)) and i G 5_. If fhere exisls j ^ 5_ such lhat Xj < Xi, fhen 
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swapping j and i would increase the magnitude of the barycenter of the group that i belongs to and so the 
objective. Now for (P_ {j,q))a feasible problem, let ,..., be its optimal partition whose corresponding 
barycenters are in ascending order and Xi be the smallest value of x in Qg, then necessarily Qi,..., Gq-i 
is optimal to solve P-{i — l,q — 1). We order therefore the values of x in ascending order and use the 
following dynamic program to compute /_, 

f-{j,q)= max f_{i-l,q-l) + {j-i + l)^i{xi,...,Xjf, (15) 

Q<i-<3 

where fj,{xi,... ,Xj) = jxjqry ^anbe computed in constant time using that 

. X Xi + {j -i)fi{xi+i,...,Xj) 

f,{xi,...,xg) = - j-i + l - 

By convention /_ {j, q) = —oo if (15) and so (P_ {j, q)) are not feasible. /„ is initialized as a grid of A: + 1 
and Q + 1 columns such that /_(0, g) = 0 for any q, f-{j, 0) = 0 and /_(j, 1) = y7r(xi,..., Xj)^ for any 
j > 1. Values of /_ are stored to compute (14). Two auxiliary variables /_ and /i_ store respectively the 
indexes of the smallest value of x in group Gg and the barycenter of the group Gq, defined by 

= argmax f-{i - l,q - 1) + {j - i + l)fi{xi,... ,Xj)'^, 
q<i<j 

fi{xi,...,Xj)<0 

i = I-{j,q). 

/_ and ;U_ are initialized by I_(j, 1) = 1 and /r_(j, 1) = /u(xi,..., xj). The same dynamic program can 
be used to compute /+, and defined similarly as /_ and /i_, by reversing the order of the values of 
X. A grid search on /(j, q, k') = f-{j, q) + f+{k' — j, Q' — q), with Q' = min(A;', Q), gives the optimal 
balance between positive and negative barycenters. A backtrack on I_ and /_|_ finally gives the best partition 
and the projection with the associated barycenters given in and 

Each dynamic program needs only to build the best partitions for the k smallest or largest partitions so 
their complexity is in 0{k^Q). The complexity of the grid search is 0{k^Q) and the complexity of the 
backtrack is 0{Q). The overall complexity of the projection is therefore 0{k‘^Q). 

5.2. Convergence. Our theoretical convergence analysis can directly be applied to this setting for a problem 
with squared loss without regularization. The feasible set is again a union of subspaces 

w € {w : w = SZv}. 

Se{0,'=,5^1=1 

However the number of largest subspaces in terms of inclusion order is smaller. They are defined by selecting 
k features among d and partitioning these k features into Q groups so that their number is N = (^){q}- 
Using classical bounds on the binomial coefficient and (11), we have for /c > 3, Q > 3, 

N < d^k^Q’^-^. 

Our analysis thus predicts that only 

n > 36 max I^Co, ^(/clogd + Qlog(/c) + {k - Q) log((5))| 

isotropic independent sub-Gaussian samples are sufficient for the projected gradient algorithm to converge. 
It produces Q + 1 cluster of features, one being a cluster of zero features, reducing dimensionality, while 
needing roughly as many samples as non-zero features. 
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6. Numerical Experiments 


We now test our methods, first on artificial dafasefs to check their robustness to noisy data. We then 
test our algorithms for feature clustering on real data extracted from movie reviews. While our approach is 
general and applies to both features and samples, we observe that our algorithms compare favorably with 
specialized algorithms for these tasks. 

6.1. Synthetic dataset. 

6.1.1. Clustering constraint on sample points. We test the robustness of our method when the information 

of the regression problem leading to the partition of the samples lies in a few features. We generate n 
data points for f = 1,. .. ,n, with Xi G d = 8, and yi G M, divided in Q = 3 clusters 

corresponding to regression tasks with weight vectors Vq. Regression labels for points Xi in cluster Qg are 
given by yi = v^xi + rjy, where py ~ AA(0, ay). We test the robustness of the algorithms to the addition of 
noisy dimensions by completing Xi with dn dimensions of noise pd ~ AA(0, Ud)- For testing the models we 
take the difference between the true label and the best prediction such that the Mean Square Error (MSE) is 
given by 

1 

Loss{y, X,W) = — V] min (yi - v^Xif. (16) 

In q=l...,Q ^ 

i=l 

The results are reported in Table 1 where the intrinsic dimension is 10 and the proportion of dimensions 
of noise dn/{d + dn) increases. On the algorithmic side, “Oracle” refers to the least-squares fit given the 
true assignments, which can be seen as the best achievable error rate, AM refers to alternate minimization, 
PG refers to projected gradient with squared loss, CG refers to conditional gradient and RC to regression 
clustering as proposed by Zhang [2003], implemented using the Harmonic K-means formulation. PG, CG 
and RC were followed by AM refinement. 1000 points were used for training, 100 for testing. The regular¬ 
ization parameters were 5-fold cross-validated using a logarithmic grid. Noise on labels is Oy = 10“^ and 
noise on added dimensions is ad = 1. Results were averaged over 50 experiments with figures after the ± 
sign corresponding to one standard deviation. 



p = 0 

p = 0.25 

p = 0.5 

p = 0.75 

p = 0.9 

p = 0.95 

Oracle 

0.52±0.08 

0.55±0.07 

0.55±0.10 

0.58±0.09 

0.71 ±0.11 

1.17±0.18 

AM 

0.52±0.08 

0.55±0.07 

5.57±4.11 

6.93±14.39 

101.08±55.49 

133.48±52.20 

PG 

1.53±7.13 

3.98±17.65 

3.20± 13.23 

5.64±20.50 

91.33±39.32 

131.48±50.90 

CG 

0.87±2.45 

1.16±4.29 

3.64± 11.02 

5.43±14.33 

91.19±53.00 

136.57±58.60 

RC 

0.52±0.08 

0.55±0.07 

5.59±20.27 

13.45±28.76 

59.19±37.97 

135.77±66.96 


Table 1. Test MSE given by (16) along proportion of added dimensions of noise p = 

dn/{d + dn)- 


All algorithms perform similarly, RC and AM get better results without added noise. None of the present 
algorithms get a significantly better behavior with a majority of noisy dimensions. 

6.1.2. Clustering constraint on features. We test the robustness of our method when with the number of 
training samples or the level of noise in the labels. We generate n data points {xi, yi) for i = 1,... ,n with 
Xi G M'^, d = 100, and y^ G M. Regression weights w have only 5 different values Ug for q = 1,..., 5, 
uniformly distributed around 0. Regression labels are given by yi = uP" Xj +p, where p rsj A/'(0,cj 2). We 
vary the number of samples n or the level of noise a and measure \\w^ — rh|| 2 > the I 2 norm of the difference 
between the true vector of weights w* and the estimated ones w. 

In Table 2 and 3, we compare the proposed algorithms to Feast Squares (ES), Feast Squares followed by 
K-means on the weights (using associated centroids as predictors) (ESK) and OSCAR [Bondell and Reich, 
2008]. For OSCAR we used a submodular approach [Bach et ah, 2012] to compute the corresponding 
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proximal algorithm, which makes it scalable. “Oracle” refers to the Least Square solution given the true 
assignments of features and can be seen as the best achievable error rate. Here too, PG refers to projected 
gradient with squared loss (initialized with the solution of Least Square followed by k-means), CG refers 
to conditional gradient, CGPG refers to conditional gradient followed by PG. When varying the number 
of samples, noise on labels is set to cr = 0.5 and when varying level of noise a number of samples is set 
to re = 150. Parameters of the algorithms were all cross-validated using a logarithmic grid. Results were 
averaged over 50 experiments and figures after the ± sign correspond to one standard deviation. 



re = 50 

re = 75 

re = 100 

re= 125 

re= 150 

Oracle 

0.16±0.06 

0.14±0.04 

0.10±0.04 

0.10±0.04 

0.09±0.03 

LS 

6L94±17.63 

5L94±16.01 

21.41 ±9.40 

L02±0.18 

0.70±0.09 

LSK 

62.93±18.05 

57.78±17.03 

10.18±14.96 

0.31±0.19 

0.19±0.12 

PG 

63.31±18.24 

52.72±16.51 

5.52±14.33 

0.14±0.09 

0.09±0.04 

CG 

6L81±17.78 

52.59±16.58 

17.24±13.87 

L20±L38 

L05±L37 

CGPG 

62.29±18.15 

50.15±17.43 

0.64±2.03 

0.15±0.19 

0.17±0.53 

OS 

61.54± 17.59 

52.87±15.90 

1L32±7.03 

L25±0.28 

0.71±0.10 


Table 2 . Measure of Hm* — w\\2, the I2 norm of the difference between the true vector of 
weights w* and the estimated ones w along number of samples re. 



cj = 0.05 

a = 0.1 

(7 = 0.5 

(7=1 

Oracle 

0.86±0.27 

L72±0.54 

8.62±2.70 

17.19±5.43 

LS 

7.04±0.92 

14.05±L82 

70.39±9.20 

140.41±18.20 

LSK 

L44±0.46 

2.88±0.91 

19.10±12.13 

48.09±27.46 

PG 

0.87±0.27 

1.74±0.52 

9.11±4.00 

26.23±18.00 

CG 

23.91±36.51 

122.31±145.77 

105.45±136.79 

155.98±177.69 

CGPG 

L52±3.13 

140.83±710.32 

17.34±53.31 

24.80± 16.32 

OS 

14.43±2.45 

18.89±3.46 

7L00±10.12 

140.33±18.83 


Table 3 . Measure of — rh||2, the I2 norm of the difference between the true vector of 
weights w* and the estimated ones w along level of noise a. 


We observe that both PG and CGPG give significantly better results than other methods and even reach 
the performance of the Oracle for re > d and for small a, while for re < d results are in the same range. 

6.2. Real data. 

6.2.1. Predicting ratings from reviews using groups of words. We perform “sentiment” analysis of newspa¬ 
per movie reviews. We use the publicly available dataset introduced by Pang and Lee [2005] which contains 
movie reviews paired with star ratings. We treat it as a regression problem, taking responses for y in (0,1) 
and word frequencies as covariates. The corpus contains re = 5006 documents and we reduced the initial 
vocabulary to d = 5623 words by eliminating stop words, rare words and words with small TF-IDF mean on 
whole corpus. We evaluate our algorithms for regression with clustered features against standard regression 
approaches: Least-Squares (LS), and Least-Squares followed by k-means on predictors (LSK), Lasso and 
Iterative Hard Thresholding (IHT). We also tested our projected gradient with sparsity constraint, initialized 
by the solution of LSK (PGS) or by the solution of CG (CGPGS). Number of clusters, sparsity constraints 
and regularization parameters were 5-fold cross-validated using respectively grids going from 5 to 15, d/2 
to d/5 and logarithmic grids. Cross validation and training were made on 80% on the dataset and tested 
on the remaining 20% it gave Q = 15 number of clusters and d/2 sparsity constraint for our algorithms. 
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Results are reported in Table 4, figures after the ± sign correspond to one standard deviation when varying 
the training and test sets on 20 experiments. 

All methods perform similarly except IHT and Lasso whose hypotheses does not seem appropriate for the 
problem. Our approaches have the benefit to reduce dimensionality from 5623 to 15 and provide meaningful 
cluster of words. The clusters with highest absolute weights are also the ones with smallest number of 
words, which confirms fhe infuifion fhaf only a few words are very discriminative. We illustrate this in 
Table 5, picking randomly words of the four clusters within which associated predictor weights Vq have 
largest magnitude. 


ES 

ESK 

PG 

CG 

CGPG 

OS 

1.51 ±0.06 

1.53±0.06 

1.52±0.06 

1.58±0.07 

1.49±0.08 

1.47 ±0.07 


PGS 

CGPGS 

IHT 

Easso 


1.53±0.06 

1.49±0.07 

2.19±0.12 

3.77±0.17 


Table 4. 


100 X mean square errors for predicting movie ratings associated with reviews. 


First and Second Cluster 

(negative) 
sizes 1 and 7 

bad, awful, 

worst, boring, ridiculous, 
watchable, suppose, disgusting. 

Last and Before Last Cluster 

(positive) 
sizes 4 and 40 

perfect,hilarious,fascinating,great 
wonderfully,perfectly,goodspirited, 
world, intelligent,wonderfully,unexpected,gem,recommendation, 
excellent,rare,unique,marvelous,good-spirited, 
mature,send,delightful,funniest 


Table 5. Clustering of words on movie reviews. We show clusters of words within which 
associated predictor weights Vq have largest magnitude. First and second one are associated 
to a negative coefficient and therefore bad feelings about movies, last and before last ones 
to a positive coefficient and good feelings about movies. 
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7. Appendix 


7.1. Formulations for classification . We present here formulations of clustering either features or samples 
when our task is to classify samples into K classes. For both settings we assume that n sample points are 
given, represented by the matrix X = (xi,..., G and corresponding labels Y = {yi,... ,yK) G 

{0,ir^-^. 

7.1.1. Clustering features for classification. Here we search K predictors W = {wi, ..., wk), each of 
them having features clustered in Q groups {Qi ,..., Qq} such that for any k, wi = vl if feature j is in 
group q. Partition of the features is shared by all predictors but each has different centroids represented in 
the vector Vk- Using an assignment matrix Z and the matrix of centroids V = {vi,, Vk), our problem can 
therefore be written 


minimize ^ lU) + |||1U||| 

subject to W = ZV, Z e {0, Z1 = 1 

in variables W G V G and Z. loss [yi,xfW) is a squared or logistic multiclass loss and 

regularization can either be seen as a standard £2 regularization on the Wk or a weighted regularization on 
the centroids Vk. 

7.1.2. Clustering samples for classification. Here our objective is to form Q groups {Qi ,..., Qq} of sam¬ 
ple points to maximize the within-group prediction performance. For classification, within each group Qq, 
samples are predicted using a common matrix of predictors = {vl ,..., u^). Our problem can be written 

1 A 

minimize - ^ loss {yi,xjv'^) + - '^Sq\\V'^\\p (17) 

i&Gq q=t 

in the variables V = {V^,..., U®) G g _ ^ such that ^ is a partition of the n 

samples, loss {yi, xfW) is a squared or logistic multiclass and | is a weighted regulariza¬ 

tion. Using an assignment matrix Z G {0, l}"x<3 ^nd auxiliary variables {W^,..., lU") G 
that VF® = if i G Qq, problem (17) can be rewritten 

minimize ^ YJi=i loss {yi, lF®^Xi^ -h f YJi=i W^\?f 
subject to W'^ = ZV^, Z G {0, l}«xQ^ Z1 = 1, 

in the variables IF G M-^xKxn^ y ^ j^dxKxQ ^ ^ (Vec(lFi),..., Vec(lF’®)), V = 

(Vec(F^),..., Vec(F‘^)) and for a matrix A, Vec(A) concatenates its columns into one vector. 

Remark that in that case we must have K > Q otherwise we output more possible answers than classes 
(in that case the problem is ill-posed). 

7.2. Clustered multitask . Our framework applies also to transfer learning by clustering similar tasks. 
Given a set of K supervised tasks like regression or binary classification, transfer learning aims at jointly 
solving these tasks, hoping that each task can benefit from the information given by other tasks. For sim¬ 
plicity, we illustrate the case of multi-category classification, which can be extended to the general multitask 
setting. When performing classification with one-versus-all majority vote, we train one binary classifier for 
each class vs. all ofhers. Using a regularizing penally such as Ihe squared £2 norm, Ihe problem of mullilask 
learning can be casl as 

minimize ^ Ef=i ELi loss(yf, rn^x*) + A (18) 

in fhe mafrix variable IF = {wi,... ,Wk) G of classifier vectors (one column per lask). We wrife 

Loss(y, X, IF) and i?(lF) fhe firsf and second term of fhis problem. Various slrafegies are used to leverage 
Ihe informalion coming from relaled fasks, such as low rank Argyriou ef al. [2008] or slrucfured norm 
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penalties Ciliberto et al. [2015] on the matrix of classifiers W. Here we follow the clustered multitask 
setting introduced in Jacob et al. [2009]. Namely we add a penalty on the classifiers (wi ,..., wk) which 
enforce them to be clustered in Q groups Qi,, Qq around centroids V = (vi,. .. ,vq) € This 

penalty can be decomposed in 

• A measure of the norm of the barycenter of centers v = ^ ^q'^q 

^mean{V) = ^KM\l 

• A measure of the variance between clusters 

Afe ^ 

^betweeni^) — ^ ^ '^(jl ^ 1 12 

• A measure of the variance within clusters 

A ® 

^u,itMn{W,V) = ^ 

9 = 1 ieGq 

The total penalty Q.{W, V) = ^mean{V) + ^between{V) + ^within{W, V) is illustrated in Figure 2. 


o 


mean 



Figure 2. Decomposed clustering penalty on K classes in the space of classifier vectors. 


The clustered multitask learning problem can then be written using an assignment matrix Z and an auxil¬ 
iary variable W Denoting If = I — the centering matrix of the K classes, we develop each term of the 
penalty, 

n^ean{V,Z) = ^ Tr{VZ^{1 - li)ZV^), 

VLbetn,een{V,Z) = ^ Tr{VZ^IiZV^), 

n^ithiniW,V,Z) = ^\\w-vz'^\\l. 

Using W = VZ^ the total penalty can then be written 

D(W,1U) = ^Tr(lU(I-n)lU'^) + ^Tr(lUniU^) + ^||lU-lU|||, 
and the problem is 

minimize Loss(y, X, W) + R{W) + D(fU, W) 

s.t.W'^ = ZV'^, Zg{0, Z1 = 1, 

in variables W £ R^^^, V gR^^^ and Z. 

19 


7.3. Convex relaxations formulations. 


7.3.1. Clustering samples for regression task. We use a squared loss Z(y, y) = 2 ( 2 /—in (3) and minimize 
in V to get a clustering problem that we can tackle using Frank-Wolfe method. We fix a partition Q and 
define for each group Cg = {ki ,..., C {1,..., d}, the matrix E G {0, that picks the Sg points 

of Qq, i.e. {Eq)ij = \ if j = ki and 0 otherwise. Therefore pg = EgU G and Xq = EqX G 
are respectively the vector of labels and the matrix of sample vectors of the group Qq. We naturally have 
EqEq = I as rows of Eq are orthonormal and Eq Eg is a diagonal matrix where Zg = diag(i7^i7g) G 
{0,1}” is the assignment vector in group Qq, i.e. {Zq)j = 1 if j ^ Qq and 0 otherwise. Z = (Zi,..., Zq) 
is therefore an assignment matrix for the partition Q. 

Minimizing in v and using the Sherman-Woodbury-Morrison formula, we obtain a function of the parti¬ 
tion 


1 Q ^ Q 


9=1 


9=1 


1 

2 n 

1 

2 n 


Q 

E - ygMs,><nl + XjXq)-^X^yq 
9=1 


9=1 


Formulating terms of the sum as solutions of an optimization problem, we get 

g=l 

1 1 

= 7^ max ^-XqX^)aq + 2yJaq, 

2n a={ai;...-,aQ) ^ ^ SqXn ^ ^ 

a,eM'’9 

where (ai;...; aq) = (af,..., a'q)'^ stacks vectors ag in one vector of size Using that 

E = {Ei ;...; Eg) = {Ef,..., Eq)^ G {0,1}’^^’^ is an orthonormal matrix, we make the change of 
variable /3 = E'^a (and so a = S/3) such that for a = (ai;...; ag), ag G M®'’, ag = Eql3. Decomposing 
Xq and pq and using E'^Eg = diag(Zq), we get 

1 1 

g=l 1 

1 1 

= — max V -d^Ej (I + -—EqXX^ET)Eql3 + 2y^SjS„/3 

2n /3eM" ^ ^ ^ SgXn ^ g J qy y g qy 

g=l 

1 Q , 

= — max Y]-Y diag(Zq)/3-diag(Zg)XX^ diag(Zq)/3 -f 2y^ diag(Zg)/3. 

2n /3elR" SgXn 

g=i 'J 

For q fixed, ( — diag(Zq)X3f^ diag(Zq)) = —xfxj if {fj) £ Oq and 0 otherwise. So 

\ *9 / ij <1 

Q 

Y — diag{Zq)XX^ diag(Zg) = XX'^ o M, 

1 ^9 
q=l 
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where M = Z{Z'^Z) ^Z'^ is the normalized equivalence matrix of the partition Q and o denotes the 
Hadamard product. Using we finally get a function of the equivalence matrix 

'tp(M) = — max o M)f3 + 2y^/3 

2 n / 3 eiR" An 

Its gradient is given by 

° (<‘+° ° ■ 

Algorithm 2 can be applied to minimize The linear oracle can indeed be computed with k-means using 
that the gradient is negative semi-definite. For a fixed Z, the linear predictors Vq for each cluster of points 
are given by 

Vq = {n\Sql +XjXq)~^X'^yq 

= {nXsql + X'^E'^EqXy^X^E'^Eqy 
= {nXsql + x'^ diag{Zq)X)~^X'^ diag{Zq)y. 

7.3.2. Convex relaxations for classification. We observe that convex relaxations for classification derive 
from computations of the convex relaxations for regression by replacing vector of labels y by the corre¬ 
sponding matrix of labels Y. 
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